home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / histogram / test2d.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  17.4 KB  |  756 lines
  1. /* histogram/test2d.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <stdio.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_math.h>
  25. #include <gsl/gsl_machine.h>
  26. #include <gsl/gsl_histogram2d.h>
  27. #include <gsl/gsl_test.h>
  28. #include <gsl/gsl_ieee_utils.h>
  29.  
  30. #define M 107
  31. #define N 239
  32. #define M1 17
  33. #define N1 23
  34. #define MR 10
  35. #define NR 5
  36.  
  37. int
  38. main (void)
  39. {
  40.   double xr[MR + 1] =
  41.     { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 };
  42.  
  43.   double yr[NR + 1] = { 90.0, 91.0, 92.0, 93.0, 94.0, 95.0 };
  44.  
  45.   gsl_histogram2d *h, *h1, *g, *hr;
  46.   size_t i, j, k;
  47.  
  48.   gsl_ieee_env_setup ();
  49.  
  50.   h = gsl_histogram2d_calloc (M, N);
  51.   h1 = gsl_histogram2d_calloc (M, N);
  52.   g = gsl_histogram2d_calloc (M, N);
  53.  
  54.   gsl_test (h->xrange == 0,
  55.         "gsl_histogram2d_calloc returns valid xrange pointer");
  56.   gsl_test (h->yrange == 0,
  57.         "gsl_histogram2d_calloc returns valid yrange pointer");
  58.   gsl_test (h->bin == 0, "gsl_histogram2d_calloc returns valid bin pointer");
  59.  
  60.   gsl_test (h->nx != M, "gsl_histogram2d_calloc returns valid nx");
  61.   gsl_test (h->ny != N, "gsl_histogram2d_calloc returns valid ny");
  62.  
  63.   hr = gsl_histogram2d_calloc_range (MR, NR, xr, yr);
  64.  
  65.   gsl_test (hr->xrange == 0,
  66.         "gsl_histogram2d_calloc_range returns valid xrange pointer");
  67.   gsl_test (hr->yrange == 0,
  68.         "gsl_histogram2d_calloc_range returns valid yrange pointer");
  69.   gsl_test (hr->bin == 0,
  70.         "gsl_histogram2d_calloc_range returns valid bin pointer");
  71.  
  72.   gsl_test (hr->nx != MR, "gsl_histogram2d_calloc_range returns valid nx");
  73.   gsl_test (hr->ny != NR, "gsl_histogram2d_calloc_range returns valid ny");
  74.  
  75.   {
  76.     int status = 0;
  77.     for (i = 0; i <= MR; i++)
  78.       {
  79.     if (hr->xrange[i] != xr[i])
  80.       {
  81.         status = 1;
  82.       }
  83.       };
  84.  
  85.     gsl_test (status,
  86.           "gsl_histogram2d_calloc_range creates xrange correctly");
  87.   }
  88.  
  89.   {
  90.     int status = 0;
  91.     for (i = 0; i <= NR; i++)
  92.       {
  93.     if (hr->yrange[i] != yr[i])
  94.       {
  95.         status = 1;
  96.       }
  97.       };
  98.  
  99.     gsl_test (status,
  100.           "gsl_histogram2d_calloc_range creates yrange correctly");
  101.   }
  102.  
  103.   for (i = 0; i <= MR; i++)
  104.     {
  105.       hr->xrange[i] = 0.0;
  106.     }
  107.  
  108.   for (i = 0; i <= NR; i++)
  109.     {
  110.       hr->yrange[i] = 0.0;
  111.     }
  112.  
  113.   {
  114.     int status = gsl_histogram2d_set_ranges (hr, xr, MR + 1, yr, NR + 1);
  115.  
  116.     for (i = 0; i <= MR; i++)
  117.       {
  118.     if (hr->xrange[i] != xr[i])
  119.       {
  120.         status = 1;
  121.       }
  122.       };
  123.  
  124.     gsl_test (status, "gsl_histogram2d_set_ranges sets xrange correctly");
  125.   }
  126.  
  127.   {
  128.     int status = 0;
  129.     for (i = 0; i <= NR; i++)
  130.       {
  131.     if (hr->yrange[i] != yr[i])
  132.       {
  133.         status = 1;
  134.       }
  135.       };
  136.  
  137.     gsl_test (status, "gsl_histogram2d_set_ranges sets yrange correctly");
  138.   }
  139.  
  140.  
  141.   k = 0;
  142.   for (i = 0; i < M; i++)
  143.     {
  144.       for (j = 0; j < N; j++)
  145.     {
  146.       k++;
  147.       gsl_histogram2d_accumulate (h, (double) i, (double) j, (double) k);
  148.     };
  149.     }
  150.  
  151.   {
  152.     int status = 0;
  153.     k = 0;
  154.     for (i = 0; i < M; i++)
  155.       {
  156.     for (j = 0; j < N; j++)
  157.       {
  158.         k++;
  159.         if (h->bin[i * N + j] != (double) k)
  160.           {
  161.         status = 1;
  162.           }
  163.       }
  164.       }
  165.     gsl_test (status,
  166.           "gsl_histogram2d_accumulate writes into array correctly");
  167.   }
  168.  
  169.   {
  170.     int status = 0;
  171.     k = 0;
  172.     for (i = 0; i < M; i++)
  173.       {
  174.     for (j = 0; j < N; j++)
  175.       {
  176.         k++;
  177.         if (gsl_histogram2d_get (h, i, j) != (double) k)
  178.           status = 1;
  179.       };
  180.       }
  181.     gsl_test (status, "gsl_histogram2d_get reads from array correctly");
  182.   }
  183.  
  184.   for (i = 0; i <= M; i++)
  185.     {
  186.       h1->xrange[i] = 100.0 + i;
  187.     }
  188.  
  189.   for (i = 0; i <= N; i++)
  190.     {
  191.       h1->yrange[i] = 900.0 + i * i;
  192.     }
  193.  
  194.   gsl_histogram2d_memcpy (h1, h);
  195.  
  196.   {
  197.     int status = 0;
  198.     for (i = 0; i <= M; i++)
  199.       {
  200.     if (h1->xrange[i] != h->xrange[i])
  201.       status = 1;
  202.       };
  203.     gsl_test (status, "gsl_histogram2d_memcpy copies bin xranges correctly");
  204.   }
  205.  
  206.   {
  207.     int status = 0;
  208.     for (i = 0; i <= N; i++)
  209.       {
  210.     if (h1->yrange[i] != h->yrange[i])
  211.       status = 1;
  212.       };
  213.     gsl_test (status, "gsl_histogram2d_memcpy copies bin yranges correctly");
  214.   }
  215.  
  216.   {
  217.     int status = 0;
  218.     for (i = 0; i < M; i++)
  219.       {
  220.     for (j = 0; j < N; j++)
  221.       {
  222.         if (gsl_histogram2d_get (h1, i, j) !=
  223.         gsl_histogram2d_get (h, i, j))
  224.           status = 1;
  225.       }
  226.       }
  227.     gsl_test (status, "gsl_histogram2d_memcpy copies bin values correctly");
  228.   }
  229.  
  230.   gsl_histogram2d_free (h1);
  231.  
  232.   h1 = gsl_histogram2d_clone (h);
  233.  
  234.   {
  235.     int status = 0;
  236.     for (i = 0; i <= M; i++)
  237.       {
  238.     if (h1->xrange[i] != h->xrange[i])
  239.       status = 1;
  240.       };
  241.     gsl_test (status, "gsl_histogram2d_clone copies bin xranges correctly");
  242.   }
  243.  
  244.   {
  245.     int status = 0;
  246.     for (i = 0; i <= N; i++)
  247.       {
  248.     if (h1->yrange[i] != h->yrange[i])
  249.       status = 1;
  250.       };
  251.     gsl_test (status, "gsl_histogram2d_clone copies bin yranges correctly");
  252.   }
  253.  
  254.   {
  255.     int status = 0;
  256.     for (i = 0; i < M; i++)
  257.       {
  258.     for (j = 0; j < N; j++)
  259.       {
  260.         if (gsl_histogram2d_get (h1, i, j) !=
  261.         gsl_histogram2d_get (h, i, j))
  262.           status = 1;
  263.       }
  264.       }
  265.     gsl_test (status, "gsl_histogram2d_clone copies bin values correctly");
  266.   }
  267.  
  268.  
  269.   gsl_histogram2d_reset (h);
  270.  
  271.   {
  272.     int status = 0;
  273.  
  274.     for (i = 0; i < M * N; i++)
  275.       {
  276.     if (h->bin[i] != 0)
  277.       status = 1;
  278.       }
  279.     gsl_test (status, "gsl_histogram2d_reset zeros array correctly");
  280.   }
  281.  
  282.   gsl_histogram2d_free (h);
  283.   h = gsl_histogram2d_calloc (M1, N1);
  284.  
  285.   {
  286.  
  287.     int status = 0;
  288.  
  289.     for (i = 0; i < M1; i++)
  290.       {
  291.     for (j = 0; j < N1; j++)
  292.       {
  293.         gsl_histogram2d_increment (h, (double) i, (double) j);
  294.  
  295.         for (k = 0; k <= i * N1 + j; k++)
  296.           {
  297.         if (h->bin[k] != 1)
  298.           {
  299.             status = 1;
  300.           }
  301.           }
  302.  
  303.         for (k = i * N1 + j + 1; k < M1 * N1; k++)
  304.           {
  305.         if (h->bin[k] != 0)
  306.           {
  307.             status = 1;
  308.           }
  309.           }
  310.       }
  311.       }
  312.     gsl_test (status, "gsl_histogram2d_increment works correctly");
  313.   }
  314.  
  315.   gsl_histogram2d_free (h);
  316.   h = gsl_histogram2d_calloc (M, N);
  317.  
  318.   {
  319.     int status = 0;
  320.     for (i = 0; i < M; i++)
  321.       {
  322.     double x0 = 0, x1 = 0;
  323.     gsl_histogram2d_get_xrange (h, i, &x0, &x1);
  324.  
  325.     if (x0 != i || x1 != i + 1)
  326.       {
  327.         status = 1;
  328.       }
  329.       }
  330.     gsl_test (status,
  331.           "gsl_histogram2d_get_xlowerlimit and xupperlimit works correctly");
  332.   }
  333.  
  334.  
  335.   {
  336.     int status = 0;
  337.     for (i = 0; i < N; i++)
  338.       {
  339.     double y0 = 0, y1 = 0;
  340.     gsl_histogram2d_get_yrange (h, i, &y0, &y1);
  341.  
  342.     if (y0 != i || y1 != i + 1)
  343.       {
  344.         status = 1;
  345.       }
  346.       }
  347.     gsl_test (status,
  348.           "gsl_histogram2d_get_ylowerlimit and yupperlimit works correctly");
  349.   }
  350.  
  351.  
  352.   {
  353.     int status = 0;
  354.     if (gsl_histogram2d_xmax (h) != M)
  355.       status = 1;
  356.     gsl_test (status, "gsl_histogram2d_xmax works correctly");
  357.   }
  358.  
  359.   {
  360.     int status = 0;
  361.     if (gsl_histogram2d_xmin (h) != 0)
  362.       status = 1;
  363.     gsl_test (status, "gsl_histogram2d_xmin works correctly");
  364.   }
  365.  
  366.   {
  367.     int status = 0;
  368.     if (gsl_histogram2d_nx (h) != M)
  369.       status = 1;
  370.     gsl_test (status, "gsl_histogram2d_nx works correctly");
  371.   }
  372.  
  373.   {
  374.     int status = 0;
  375.     if (gsl_histogram2d_ymax (h) != N)
  376.       status = 1;
  377.     gsl_test (status, "gsl_histogram2d_ymax works correctly");
  378.   }
  379.  
  380.   {
  381.     int status = 0;
  382.     if (gsl_histogram2d_ymin (h) != 0)
  383.       status = 1;
  384.     gsl_test (status, "gsl_histogram2d_ymin works correctly");
  385.   }
  386.  
  387.   {
  388.     int status = 0;
  389.     if (gsl_histogram2d_ny (h) != N)
  390.       status = 1;
  391.     gsl_test (status, "gsl_histogram2d_ny works correctly");
  392.   }
  393.  
  394.   h->bin[3 * N + 2] = 123456.0;
  395.   h->bin[4 * N + 3] = -654321;
  396.  
  397.   {
  398.     double max = gsl_histogram2d_max_val (h);
  399.     gsl_test (max != 123456.0, "gsl_histogram2d_max_val finds maximum value");
  400.   }
  401.  
  402.   {
  403.     double min = gsl_histogram2d_min_val (h);
  404.     gsl_test (min != -654321.0,
  405.           "gsl_histogram2d_min_val finds minimum value");
  406.   }
  407.  
  408.   {
  409.     size_t imax, jmax;
  410.     gsl_histogram2d_max_bin (h, &imax, &jmax);
  411.     gsl_test (imax != 3
  412.           || jmax != 2,
  413.           "gsl_histogram2d_max_bin finds maximum value bin");
  414.   }
  415.  
  416.   {
  417.     size_t imin, jmin;
  418.     gsl_histogram2d_min_bin (h, &imin, &jmin);
  419.     gsl_test (imin != 4
  420.           || jmin != 3, "gsl_histogram2d_min_bin find minimum value bin");
  421.   }
  422.  
  423.   for (i = 0; i < M * N; i++)
  424.     {
  425.       h->bin[i] = i + 27;
  426.       g->bin[i] = (i + 27) * (i + 1);
  427.     }
  428.  
  429.   {
  430.     double sum = gsl_histogram2d_sum (h);
  431.     gsl_test (sum != N * M * 27 + ((N * M - 1) * N * M) / 2,
  432.           "gsl_histogram2d_sum sums all bin values correctly");
  433.   }
  434.  
  435.   {
  436.     /* first test... */
  437.     const double xpos = 0.6;
  438.     const double ypos = 0.85;
  439.     double xmean;
  440.     double ymean;
  441.     size_t xbin;
  442.     size_t ybin;
  443.     gsl_histogram2d *h3 = gsl_histogram2d_alloc (M, N);
  444.     gsl_histogram2d_set_ranges_uniform (h3, 0, 1, 0, 1);
  445.     gsl_histogram2d_increment (h3, xpos, ypos);
  446.     gsl_histogram2d_find (h3, xpos, ypos, &xbin, &ybin);
  447.     xmean = gsl_histogram2d_xmean (h3);
  448.     ymean = gsl_histogram2d_ymean (h3);
  449.  
  450.     {
  451.       double expected_xmean = (h3->xrange[xbin] + h3->xrange[xbin + 1]) / 2.0;
  452.       double expected_ymean = (h3->yrange[ybin] + h3->yrange[ybin + 1]) / 2.0;
  453.       gsl_test_abs (xmean, expected_xmean, 100.0 * GSL_DBL_EPSILON,
  454.             "gsl_histogram2d_xmean works correctly");
  455.       gsl_test_abs (ymean, expected_ymean, 100.0 * GSL_DBL_EPSILON,
  456.             "gsl_histogram2d_ymean works correctly");
  457.     };
  458.     gsl_histogram2d_free (h3);
  459.   }
  460.  
  461.   {
  462.     /* test it with bivariate normal distribution */
  463.     const double xmean = 0.7;
  464.     const double ymean = 0.7;
  465.     const double xsigma = 0.1;
  466.     const double ysigma = 0.1;
  467.     const double correl = 0.5;
  468.     const double norm =
  469.       10.0 / M_PI / xsigma / ysigma / sqrt (1.0 - correl * correl);
  470.     size_t xbin;
  471.     size_t ybin;
  472.     gsl_histogram2d *h3 = gsl_histogram2d_alloc (M, N);
  473.     gsl_histogram2d_set_ranges_uniform (h3, 0, 1, 0, 1);
  474.     /* initialize with 2d gauss pdf in two directions */
  475.     for (xbin = 0; xbin < M; xbin++)
  476.       {
  477.     double xi =
  478.       ((h3->xrange[xbin] + h3->xrange[xbin + 1]) / 2.0 - xmean) / xsigma;
  479.     for (ybin = 0; ybin < N; ybin++)
  480.       {
  481.         double yi =
  482.           ((h3->yrange[ybin] + h3->yrange[ybin + 1]) / 2.0 -
  483.            ymean) / ysigma;
  484.         double prob =
  485.           norm * exp (-(xi * xi - 2.0 * correl * xi * yi + yi * yi) /
  486.               2.0 / (1 - correl * correl));
  487.         h3->bin[xbin * N + ybin] = prob;
  488.       }
  489.       }
  490.     {
  491.       double xs = gsl_histogram2d_xsigma (h3);
  492.       double ys = gsl_histogram2d_ysigma (h3);
  493.       /* evaluate results and compare with parameters */
  494.  
  495.       gsl_test_abs (gsl_histogram2d_xmean (h3), xmean, 2.0/M,
  496.                     "gsl_histogram2d_xmean works correctly");
  497.       gsl_test_abs (gsl_histogram2d_ymean (h3), ymean, 2.0/N,
  498.                     "gsl_histogram2d_ymean works correctly");
  499.       gsl_test_abs (xs, xsigma, 2.0/M,
  500.                     "gsl_histogram2d_xsigma works correctly");
  501.       gsl_test_abs (ys, ysigma, 2.0/N,
  502.                     "gsl_histogram2d_ysigma works correctly");
  503.       gsl_test_abs (gsl_histogram2d_cov (h3) / xs / ys, correl,
  504.                     2.0/((M < N) ? M : N),
  505.                     "gsl_histogram2d_cov works correctly");
  506.     }
  507.     gsl_histogram2d_free (h3);
  508.   }
  509.  
  510.   gsl_histogram2d_memcpy (h1, g);
  511.   gsl_histogram2d_add (h1, h);
  512.  
  513.   {
  514.     int status = 0;
  515.     for (i = 0; i < M * N; i++)
  516.       {
  517.     if (h1->bin[i] != g->bin[i] + h->bin[i])
  518.       status = 1;
  519.       }
  520.     gsl_test (status, "gsl_histogram2d_add works correctly");
  521.   }
  522.  
  523.   gsl_histogram2d_memcpy (h1, g);
  524.   gsl_histogram2d_sub (h1, h);
  525.  
  526.   {
  527.     int status = 0;
  528.     for (i = 0; i < M * N; i++)
  529.       {
  530.     if (h1->bin[i] != g->bin[i] - h->bin[i])
  531.       status = 1;
  532.       }
  533.     gsl_test (status, "gsl_histogram2d_sub works correctly");
  534.   }
  535.  
  536.  
  537.   gsl_histogram2d_memcpy (h1, g);
  538.   gsl_histogram2d_mul (h1, h);
  539.  
  540.   {
  541.     int status = 0;
  542.     for (i = 0; i < M * N; i++)
  543.       {
  544.     if (h1->bin[i] != g->bin[i] * h->bin[i])
  545.       status = 1;
  546.       }
  547.     gsl_test (status, "gsl_histogram2d_mul works correctly");
  548.   }
  549.  
  550.   gsl_histogram2d_memcpy (h1, g);
  551.   gsl_histogram2d_div (h1, h);
  552.  
  553.   {
  554.     int status = 0;
  555.     for (i = 0; i < M * N; i++)
  556.       {
  557.     if (h1->bin[i] != g->bin[i] / h->bin[i])
  558.       status = 1;
  559.       }
  560.     gsl_test (status, "gsl_histogram2d_div works correctly");
  561.   }
  562.  
  563.   gsl_histogram2d_memcpy (h1, g);
  564.   gsl_histogram2d_scale (h1, 0.5);
  565.  
  566.   {
  567.     int status = 0;
  568.     for (i = 0; i < M * N; i++)
  569.       {
  570.     if (h1->bin[i] != 0.5 * g->bin[i])
  571.       status = 1;
  572.       }
  573.     gsl_test (status, "gsl_histogram2d_scale works correctly");
  574.   }
  575.  
  576.   gsl_histogram2d_memcpy (h1, g);
  577.   gsl_histogram2d_shift (h1, 0.25);
  578.  
  579.   {
  580.     int status = 0;
  581.     for (i = 0; i < M * N; i++)
  582.       {
  583.     if (h1->bin[i] != 0.25 + g->bin[i])
  584.       status = 1;
  585.       }
  586.     gsl_test (status, "gsl_histogram2d_shift works correctly");
  587.   }
  588.  
  589.   gsl_histogram2d_free (h);    /* free whatever is in h */
  590.  
  591.   h = gsl_histogram2d_calloc_uniform (M1, N1, 0.0, 5.0, 0.0, 5.0);
  592.  
  593.   gsl_test (h->xrange == 0,
  594.         "gsl_histogram2d_calloc_uniform returns valid range pointer");
  595.   gsl_test (h->yrange == 0,
  596.         "gsl_histogram2d_calloc_uniform returns valid range pointer");
  597.   gsl_test (h->bin == 0,
  598.         "gsl_histogram2d_calloc_uniform returns valid bin pointer");
  599.   gsl_test (h->nx != M1, "gsl_histogram2d_calloc_uniform returns valid nx");
  600.   gsl_test (h->ny != N1, "gsl_histogram2d_calloc_uniform returns valid ny");
  601.  
  602.   gsl_histogram2d_accumulate (h, 0.0, 3.01, 1.0);
  603.   gsl_histogram2d_accumulate (h, 0.1, 2.01, 2.0);
  604.   gsl_histogram2d_accumulate (h, 0.2, 1.01, 3.0);
  605.   gsl_histogram2d_accumulate (h, 0.3, 0.01, 4.0);
  606.  
  607.   {
  608.     size_t i1, i2, i3, i4;
  609.     size_t j1, j2, j3, j4;
  610.     double expected;
  611.     int status;
  612.     status = gsl_histogram2d_find (h, 0.0, 3.01, &i1, &j1);
  613.     status = gsl_histogram2d_find (h, 0.1, 2.01, &i2, &j2);
  614.     status = gsl_histogram2d_find (h, 0.2, 1.01, &i3, &j3);
  615.     status = gsl_histogram2d_find (h, 0.3, 0.01, &i4, &j4);
  616.  
  617.     for (i = 0; i < M1; i++)
  618.       {
  619.     for (j = 0; j < N1; j++)
  620.       {
  621.         if (i == i1 && j == j1)
  622.           {
  623.         expected = 1.0;
  624.           }
  625.         else if (i == i2 && j == j2)
  626.           {
  627.         expected = 2.0;
  628.           }
  629.         else if (i == i3 && j == j3)
  630.           {
  631.         expected = 3.0;
  632.           }
  633.         else if (i == i4 && j == j4)
  634.           {
  635.         expected = 4.0;
  636.           }
  637.         else
  638.           {
  639.         expected = 0.0;
  640.           }
  641.  
  642.         if (h->bin[i * N1 + j] != expected)
  643.           {
  644.         status = 1;
  645.           }
  646.       }
  647.       }
  648.     gsl_test (status, "gsl_histogram2d_find works correctly");
  649.   }
  650.  
  651.   {
  652.     FILE *f = fopen ("test.txt", "w");
  653.     gsl_histogram2d_fprintf (f, h, "%.19g", "%.19g");
  654.     fclose (f);
  655.   }
  656.  
  657.   {
  658.     FILE *f = fopen ("test.txt", "r");
  659.     gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1);
  660.     int status = 0;
  661.  
  662.     gsl_histogram2d_fscanf (f, hh);
  663.  
  664.     for (i = 0; i <= M1; i++)
  665.       {
  666.     if (h->xrange[i] != hh->xrange[i])
  667.       {
  668.         printf ("xrange[%d] : %g orig vs %g\n",
  669.             (int) i, h->xrange[i], hh->xrange[i]);
  670.         status = 1;
  671.       }
  672.       }
  673.  
  674.     for (j = 0; j <= N1; j++)
  675.       {
  676.     if (h->yrange[j] != hh->yrange[j])
  677.       {
  678.         printf ("yrange[%d] : %g orig vs %g\n",
  679.             (int) j, h->yrange[j], hh->yrange[j]);
  680.         status = 1;
  681.       }
  682.       }
  683.  
  684.     for (i = 0; i < M1 * N1; i++)
  685.       {
  686.     if (h->bin[i] != hh->bin[i])
  687.       {
  688.         printf ("bin[%d] : %g orig vs %g\n",
  689.             (int) i, h->bin[i], hh->bin[i]);
  690.         status = 1;
  691.       }
  692.       }
  693.  
  694.     gsl_test (status, "gsl_histogram2d_fprintf and fscanf work correctly");
  695.  
  696.     gsl_histogram2d_free (hh);
  697.     fclose (f);
  698.   }
  699.  
  700.   {
  701.     FILE *f = fopen ("test.dat", "wb");
  702.     gsl_histogram2d_fwrite (f, h);
  703.     fclose (f);
  704.   }
  705.  
  706.   {
  707.     FILE *f = fopen ("test.dat", "rb");
  708.     gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1);
  709.     int status = 0;
  710.  
  711.     gsl_histogram2d_fread (f, hh);
  712.  
  713.     for (i = 0; i <= M1; i++)
  714.       {
  715.     if (h->xrange[i] != hh->xrange[i])
  716.       {
  717.         printf ("xrange[%d] : %g orig vs %g\n",
  718.             (int) i, h->xrange[i], hh->xrange[i]);
  719.         status = 1;
  720.       }
  721.       }
  722.  
  723.     for (j = 0; j <= N1; j++)
  724.       {
  725.     if (h->yrange[j] != hh->yrange[j])
  726.       {
  727.         printf ("yrange[%d] : %g orig vs %g\n",
  728.             (int) j, h->yrange[j], hh->yrange[j]);
  729.         status = 1;
  730.       }
  731.       }
  732.  
  733.     for (i = 0; i < M1 * N1; i++)
  734.       {
  735.     if (h->bin[i] != hh->bin[i])
  736.       {
  737.         printf ("bin[%d] : %g orig vs %g\n",
  738.             (int) i, h->bin[i], hh->bin[i]);
  739.         status = 1;
  740.       }
  741.       }
  742.  
  743.     gsl_test (status, "gsl_histogram2d_fwrite and fread work correctly");
  744.  
  745.     gsl_histogram2d_free (hh);
  746.     fclose (f);
  747.   }
  748.  
  749.   gsl_histogram2d_free (h);
  750.   gsl_histogram2d_free (h1);
  751.   gsl_histogram2d_free (g);
  752.   gsl_histogram2d_free (hr);
  753.  
  754.   exit (gsl_test_summary ());
  755. }
  756.